*	See Gibbs estimation of microstructure models: Teaching notes
	November, 2006
	Joel Hasbrouck
____________________________________________________________________________________________________

	RollGibbs Analyze q.sas

____________________________________________________________________________________________________;

*	Modify the following statement to set current directory to location of program;
x 'cd c:\Active\Mkt Microstructure Book\Programs Fall 2006';
options nodate nocenter nonumber ps=70 ls=120  sasautos=('./Sas Macros','.') mprint; 
libname SASData '.\Datasets';
libname this '.';

/*ods pdf file='./RollGibbs 01.pdf' style=serifPrinter;*/

%let Infinity=1e30;
%let eps=1e-30;

*___________________________________________________________________________________________________

	Generate some price data
____________________________________________________________________________________________________;
%let nObs=20;
%let sdu = .01;
*	Three values of c are considered.;
%let c = .01;
%let c = 0.005;
%let c = 0.001;
data sim;
	m = 0;
	do t=1 to &nObs;
		m = m + &sdu*normal(345699);
		q = 1-2*(uniform(246)>.5);
		p = m + &c*q;
		keep q p t;
		output;
	end;
	run;
goptions ftext=complex htext=2 colors=(black);
symbol1 v=dot i=join w=1;
axis1 minor=none label=('Time');
axis2 minor=none label=(a=90 'Price');
proc gplot data=sim;
	plot p*t=1 / haxis=axis1 vaxis=axis2;
	title "Simulated price path with c=&c";
	run;
	quit;

%let nSweeps=10000;

*___________________________________________________________________________________________________

	Simulate parameters and q's
____________________________________________________________________________________________________;
%RollGibbs(sim, nSweeps=&nSweeps, cStart=&c, regDraw=0, varuStart=%sysevalf(&sdu*&sdu), varuDraw=0);



*___________________________________________________________________________________________________

	Analyze the simulated estimates
____________________________________________________________________________________________________;
%let nDrop = 0;
data qOut2;
	set qOut;
	if sweep>&nDrop then output;
	run;
proc transpose data=qOut2 out=qt prefix=q;
	var q;
	by sweep;
	run;
proc means data=qt (drop=sweep) noprint;
	output out=qtMeans (drop=_type_ _freq_) mean= ;
	run;
proc transpose data=qtMeans out=qMeans (rename=(col1=q));
	run;
data qMeans;
	set qMeans;
	t = _n_;
	keep t q;
	run;
proc print data=qMeans;
	run;
data qBoth;
	merge qMeans (rename=(q=qEst)) sim;
	by t;
	run;
symbol1 i=needle v=none w=5;
symbol2 i=none v=dot h=1.5;
axis2 order=(-1 to 1 by 1) label=(a=90 "E[q(t)]") minor=none value=("-1 (Sell)" " " "+1 (Buy)");
proc gplot data=qBoth;
	plot qEst*t=1 q*t=2 / haxis=axis1 vaxis=axis2 overlay;
	title "Actual and simulated trade directions with c=&c";
	run;
	quit;
